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ABSTRACT 

We aim at reproducing the height dependence of sunspot wave signatures obtained from spectropo- 
larimetric observations through 3D MHD numerical simulations. A magneto-static sunspot model 
based on the properties of the observed sunspot is constructed and perturbed at the photosphere 
introducing the fluctuations measured with the Si I A 10827 A line. The results of the simulations are 
compared with the oscillations observed simultaneously at different heights from the He I A 10830 A 
line, the Call H core and the Fel blends in the wings of the Call H line. The simulations show a 
remarkable agreement with the observations. They reproduce the velocity maps and power spectra 
at the formation heights of the observed lines, as well as the phase and amplification spectra between 
several pair of lines. We find that the stronger shocks at the chromosphere are accompanied with a 
delay between the observed signal and the simulated one at the corresponding height, indicating that 
shocks shift the formation height of the chromospheric lines to higher layers. Since the simulated wave 
propagation matches very well the properties of the observed one, we are able to use the numerical 
calculations to quantify the energy contribution of the magneto-acoustic waves to the chromospheric 
heating in sunspots. Our findings indicate that the energy supplied by these waves is too low to 
balance the chromospheric radiative losses. The energy contained at the formation height of the low- 
ermost Si I A 10827 A line in the form of slow magneto-acoustic waves is already insufficient to heat 
the higher layers, and the acoustic energy which reaches the chromosphere is around 3-9 times lower 
than the required amount of energy. The contribution of the magnetic energy is even lower. 
Subject headings: MHD; Sun: chromosphere; Sun: oscillations; Sun: photosphere; sunspots 



1. INTRODUCTION 

The question about the processes which heat the stel- 
lar outer atmospheres is one of the most intriguing unan- 
swered problems in astrophysics. Among all the mecha- 
nisms which have been proposed to account for these en- 
ergy losses, two of them seem to be the most promising 
ones: mechanical heating by up ward-propagating waves 
generated in the convectio n zone ()Alfvenlll947l : IBiermannI 
119481 : ISchwarzschildl 11948") and Joule heating driven by 
magnetic field reconnection an d resistive dissipation o f 
electric curre nts ([P arker 1983; Hev vaerts fc Priestlll983D . 
The work by iSocas-Navarro (2005j concluded that likely 
the Joule heating mechanism cannot provide the domi- 
nant source to heat the sunspot chromosphere, leaving 
the energy transport by waves as a plausible candidate. 

The magnitude of the energy supplied to the 
quiet Solar chromosphere by high-frequency acoustic 
waves has been w i dely s tudie d over the last years. 
iFossum fc Carlssonl ()2005L l2006f ) analyzed temporal se- 
ries of the quiet Sun observations obtained with the 
Transition Region And Coronal Explorer (TRACE) in 
two continuum bands at 1600 and 1700 A. From the 
comparison of these observations with one-dimensional 
(ID) numerical simulations, they retrieved an acoustic 
flux supplied to the chromosphere by high frequency 
waves (at 5-28 mHz) of at least 10 times lower than 
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the required amount of energy to account for the chro- 
mospheric radiative losses. However, some authors have 
criticized this result . Using three-dimensional (3D ) nu- 
merical simu l ations , IWedemever-Bohm et al.l ()2007[ ) and 
iCuntz et al.l ()2007l ) have argued that the limited spa- 
tial resolution of TRACE, around 1", hides a factor 
of 10 in the short-period energy flux. Kalkofen! (j2007l ) 
shared this point of view, and claimed that the observa- 
tions of chromospheric radiation s upport the heating by 
the dissipation of acoustic waves. iCarlsson et al.l ()2007l ) 
analyzed high-resolution time sequences of Call H fil- 
tergrams from SOT/HINODE, finding a larger acous- 
tic power, but still too low to balance the chromo- 
sph eric radiative losses. A n eve n larger value was found 
by iBello Gonzalez et al.l (j2009D. Recently, using the 
IMAX spectroDolari meter (iMartinez Pillet et al.l I2OIOD 
onboar d SUNRISE (iBarthol et al I 120101 : ISolanki d!ld] 
I2OIOI ). iBello Gonzalez et al.l (|2010l ) obtained an acous- 
tic flux at a height of 250 km above the photosphere 
only around 2 times l ower than th e chrom ospheric radia- 
tive losses. However, iFleck et al.l ()2010f ) argue that en- 
ergy fluxes derived from measured high frequency waves 
may be largely overestimated, since their power may be 
caused by line formation effects in a dynamic atmosphere 
rather than propagating high frequency acoustic waves. 

With regards to quiet Sun waves with frequencies i n 
the 3 and 5 mHz bands, according to iBeck et al.l (|2009f ). 
their wave energy is insufficient to maintain the temper- 
ature increase at layers higher than 500 km above the 
photosphere. 

Some works have pointed out that the energy flux as- 
sociated to other wave modes may be added up to bal- 
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ance the chromospheric radiative losses. iStraus et al.l 
(|2QQ8f ) detected low- frequency gravity waves in the Sun's 
atmosphere, and claimed that their amount of energy 
is comparable to the radiative losses of the chromo- 
sphere. On the other hand, the simulated acoustic flux 
obtained by them is a factor of 10 smaller, consistent with 
iFossum fc Carlsso n (2005). Regarding the Alfven waves, 
Ide Pontieu et al.l (j2007.) claim to detect them in spicules 
in the upper chromosphere and estimate the energy flux 
carried by these waves to be sufficient to accelerate the 
solar wind and possibly to heat the quiet corona. 

Most of the works that address the problem of upper 
atmosphere heating by waves refer to quiet Sun regions. 
In this paper we aim to study the energy balance in a 
su nspot umbra. An ea rlier work in this topic is the one 
by iKneer et al.l (|1981f ). who analize the time lags and 
rms velocities between the lines Nal Di and Nal D2 and 
found that the waves supply to the chromosphere an en- 
ergy flux of about 5 X 10^ erg cm~^ s~^, which is much 
lower than the chromospheric radiative losses of about 
10^ erg cm~^ s~^. Our approach consists in using 3D nu- 
merical simulations to reproduce the oscillatory pattern 
from the photosphere to the chromosphere obtained from 
temporal series of spectropolarimetric data. The resem- 
blance between the simulated and observed waves allows 
us to extract conclusions from the numerical calculations 
provided that they correspond to a realistic phenomenon. 
In Section [2] we describe the observational data and the 
numerical method employed. Sections [3] and H] discuss 
the construction and properties of the magnetohydro- 
static (MHS) model of sunspot and the procedure used 
to introduce the driver, respectively. The analysis of the 
simulations and their comparison with the observations 
is performed in Section [51 while in Section [6] we use these 
results to evaluate the energy balance. Finally, the dis- 
cussion and conclusions are presented in Section [71 

2. OBSERVATIONAL AND NUMERICAL PROCEDURES 

The observational data consists of a temporal series 
of co-spatial and simultaneous data obtained with two 
different instru nients, the POlarin ietric LIttrow Spectro- 
graph (POLIS, |Beck_elay [2QQ5bh and the Tener ife In- 
frared Polarimeter II (TIP-II,[Coilados et al. 2007), both 
attached to the German Vacuum Tower telescope at the 
Observatorio del Teide at Tenerife on August 27^^ 2007. 
TIP-II provides Stokes spectra of the 10830 A region, in- 
cluding the Si I A 10827 A and Hei A 10830 A spectral 
lines, with a spatial sampling of 0''18 per pixel. The 
blue channel of POLIS yields intensity spectra of the 
Call H A 3968 A line and some photospheric Fel line 
blends in its wings, with a spatial sampling of 0''29 per 
pixel. The slit was placed over the center of a sunspot lo- 
cated near the center of the Su n. The thorough ana lysis 
of this dataset was presented in lFelipe et al.l ()2010bl ). 

We have performed 3D numerical simulations to repro- 
duce the observed wave pattern from the photosphere to 
the chromo sphere . The numeric al method is described 
in detail in [Felipe et al.l (|2010a[ ). The code solves the 
nonlinear MHD equations for perturbations. A MHS 
model of sunspot is perturbed with a driver force in 
the momentum equation. The properties of the MHS 
model and the driver used in the current work, are dis- 
cussed in the following sections. A p erfectly matched 
layer (PML) boundary condition ()Berenger..l994) is ap- 



plied to all boundaries in order to avoid wave reflections. 
The radiative energy losses are implemented following 
Newton's cooling law: 

T 

Qrad 1 (1) 

where Ti is the perturbation in the temperature, tr is 
the radiative relaxation time, and Cy is the specific heat 
at constant volume. Following iSpiegel (1957^), tr is given 

by 

= 16^^' 

where x is the mean absortion coefficient and ctr is the 
Stefan-Boltzmann constant. This expression is valid at 
photospheric heights, but not at the chromosphere as it 
was derived by Spiegel in the approximation of local ther- 
modynamic equilibrium. As the values of tr given by the 
Spiegel formula are not correct at chromospheric heights, 
we took the freedom to modify them in order to mimic 
the low value of r^^ = 10 s obtained from the observa- 
tional analysis of the w ave prop agation at chromospheric 
heights in .Feli pe et all (|2010b[ i. 

3. MHS MODEL OF THE SUNSPOT 

The MHS model is constructe d following the method 
bv iKhomenko fc Colladosl (|2008D . These authors devel- 
oped a method to calculate a thick sunspot structure 
in magnetostatic equilibrium with distributed currents, 
i.e., showing a continuous variation of field strength and 
gas pressure across the spot, from sub-photospheric to 
chromospheric layers. In current-distributed models, the 
field decreases from its value at the axis of the sunspot 
to almost zero at large radial distances. These models 
are constructed by combining the advantages of two dif- 
ferent methods: in the deep layers, the magnetic topol- 
ogy is set and the thermod ynamic variables are forced to 
match with this structure (jSchliiter fc Temesvarvl[T958l : 
lLowlll975l . ll980l ). constructing the so-called "self-similar" 
models; at photospheric heights and above, the pressure 
distribution is prescribed as boundary condition at the 
axis of the sunspot and in the distant quiet Sun atmo- 
sphere and in between the pressure and magnetic field 
are iter at ively changed until the system reaches an equi- 
librium state (|Pizzolll986l ). 

The stratification of thermodynamic variables in the 
photosphere (needed for the boundary conditions of 
the above methods) was retrieved from the inversion 
of the Si I Stokes spectra in the umbra and the quiet 
Sun atmosphere, averaged in time an d space, using SIR 
(|Ruiz Cobo del Toro Iniestd 119921 ). From the mver- 
sion of the Si I Stokes profiles we can retrieve the strat- 
ification of magnetic field and thermodynamic variables 
up to a height of almost 800 km in the atmosphere of the 
observed sunspot and its surroundings. Since the numer- 
ical wave calculations need the distribution of gas pres- 
sure deeper below the photosphere as well as at higher 
chromospheric layers, we have smoothly joined our es- 
timation of the photospheric variables with other mod- 
els from t he literature. As field-free atmosphe re we used 
model S of lChristensen-Dalsgaard et al.l () 19961 ) at deeper 
layers and the VAL-C model of the solar chromosphere 
([Vernazza et al.iri98l[ ) in t he uppe r layer s. At the axis 
of the sunspot we use the lAvrettl (|1981l ) model in the 
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upper layer s, while the deep layers w ere extracted from 
a model by iKosovichev et all (j2QQQf ) obtained from he- 
lioseismic inversions of the phase speed in sunspots. The 
resulting sunspot atmosphere was shifted down 350 km 
in the vertical direction in order to account for the esti- 
mated Wilson depression (according to horizontal pres- 
sure balance). These stratifications were set as boundary 
conditions at the quiet Sun and at the umbral axis at- 
mosphere, respectively. 

The method also requires to characterize the horizon- 
tal variations of the magnetic field of the structure at 
some height. At each spatial position of our observations 
(along the slit) we have averaged the Stokes profiles of 
the Si I A 10827 A line in time, and we have inverted the 
resulting Stokes vectors with SIR. From a gaussian fit to 
the spatial variation of the magnetic field obtained from 
the inversion at the formation height of the silicon line 
we have obtained the parameters which define the radial 
distribution of the magnetic field, i.e., its strength at the 
axis and the effective diameter. 
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Fig. 1. — MHS sunspot model constructed to match the observed 
properties. Top: Magnetic field strength; bottom: Temperature. 
White thin lines are magnetic field lines. White thick lines with 
labels are the contours of The yellow line is the layer with 

optical depth unity. The height z = Mm corresponds to the 
height where optical depth is unity at quiet Sun atmosphere. 



Figure [T] shows the distribution of the magnetic field 
and temperature in the complete sunspot model. Most 
of the magnetic field is concentrated around the axis of 
the model, inside a radius of 10 Mm, and it is weaker at 
larger distances from the center. The orientation of the 
magnetic field lines changes very fast from being verti- 
cal at the axis to almost horizontal at radial distances 
larger than 30 Mm below 2; = km. Above this layer 
the magnetic field lines spread increasing their inclina- 
tion with the distance to the axis. The properties of the 
model are shown in more detail in Figure [21 where the 
distribution with radius at several heights (left panels) 
and the stratification with height at different distances 
from the axis of the sunspot (right panels) are plotted. 
Note that we only show values corresponding to the green 
box in Figure [U as this box is used in the simulations. 
At the axis of the sunspot the magnetic field drops from 
1300 G at z = -1 Mm to 1.1 kG at z = 1 Mm. It is 
vertical at r = Mm and its inclination at photospheric 
layers increases with the radius, reaching an inclination 



around 20^ at the rightmost point of the simulation do- 
main. The gas pressure has a deficit at the magnetized 
regions around the center of the sunspot. Below z = 
km this deficit decreases with depth, as can be seen from 
the radial distribution of gas pressure at z = — 1 Mm, 
and it almost disappears at about —2 Mm depth. The 
gas pressure changes by 6 orders of magnitude from its 
value at z = — 1 Mm to z = 1 Mm. The squared ratio 
between the sound speed and the Alfven speed has very 
strong variations, from 10^ at z = —1 Mm to 10~^ at 
z = \ Mm (taking the values at the axis). 




RADIUS [Mm] 



Z [Mm] 



Fig. 2. — Distribution with radial distance (left panels) and with 
depth (right panels) of the magnetic field strength, pressure, ratio 
<^s/'^A' magnetic field inclination for the obtained sunspot 

model. The radial pressure distributions are normalized to their 
values at the right boundary (non-magnetic). Only the values cor- 
responding to heights/distances in the green box from Figure ^ are 
shown. 

It must be noted that the observed photospheric mag- 
netic field strength was around 2000 G, while the model 
has a value between 1100 G and 1300 G. This mismatch 
is the result of the complex merging of all the different 
model atmospheres. As explained above, the stratifica- 
tion of the thermodynamic variables in the sunspot axis 
has been obtained by merging three independent models: 
the one by Kosovichev et al. (subphotospheric layers), 
the one retrieved from the inversions (photospheric lay- 
ers) and Avrett's model (upper layers). Also, the pho- 
tospheric model is formed by the model of Christensen- 
Dalsgaard (subphotospheric layers), the one from the in- 
versions (photosphere) and VAL-C (upper layers). The 
different models did not coincide at the same height for 
all parameters, or even, some times, discontinuities ap- 
peared in the stratification of a given parameter, which 
made necessary an interpolation at some layers to merge 
them smoothly. At the end, the resulting stratifications 
at the sunspot axis and field-free atmospheres did not 
satisfy hydrostatic equilibrium and pressures were re- 
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calculated to impose it. This re-calculation gave rise 
to a variation of the pressure stratification at both the 
sunspot axis and quiet sun atmospheres, with a reduced 
pressure deficit (which turned at the end into a lower 
sunspot magnetic field strength). The later merging with 
Low's model to generate a deep sunspot also suffers from 
the same problem, but this represented a minor correc- 
tion. At this point, we could have opted to improve 
the sunspot model to match the observed magnetic field 
strength. But, given that the retrieved stratifications of 
the parameters which mainly determine the properties 
of wave propagation (Alfven and sound velocities and 
pressure scale height after the iterative process) are very 
similar to those obtained from the inversions (see Figure 
[3]), we were confident that the model was adequate to 
study the propagation of waves in the atmosphere of the 
observed sunspot. 
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Fig. 3. — Comparison between the computed model of sunspot 
(black solid line) and the inversion (red dashed line) at the axis of 
the sunspot for a range of heights around the photosphere. From 
top to bottom: sound velocity, Alfven velocity and pressure scale 
height. 



4. INTRODUCTION OF THE DRIVER 

In this work our aim is the reproduction of the observed 
wave pattern by means of numerical calculations. In this 
respect it is crucial to choose the most appropriate way to 
introduce the observed velocity as a driver. We have cho- 
sen the velocity measured with the Si I line as the driver 
of the simulation, since it is the line which is formed the 
deepest over the set of lines that we have observed. At 
the formation height of the Si I line, the numerical simu- 
lation should have a vertical velocity as close as possible 
to the measured LOS velocity. The photospheric oscil- 
lations are dominated by waves in the 5 minute band 
and, thus, the power excited in the simulation in this 
band must resemble the observed one. However, it is 
even more critical to introduce correctly the power at 
higher frequencies. Waves with frequencies above the 
cutoff propagate upward and dominate the higher layers. 
The wave pattern at the chromosphere will depend on the 



power introduced by the driver at those high frequencies 
as well as on the initial phase of these high frequency 
waves. 

It is interesting to note that the simulations by 
iCarlsson fc Steiril ( 19971 ) showed that it is possible to 
derive a transfer function for accurately relating the 
observed velocity with a piston velocity at the bot- 
tom boundary, where the latter may be located at a 
deeper position. This a pproach is not valid in our case. 
^Carls son fc Steiril (|1997[ ) performed one dimensional hy- 
drodynamic simulations, where they only can propagate 
acoustic longitudinal waves in the vertical direction, and 
there is a unique transfer function for this wave, without 
ambiguity. However, in the three dimensional case, there 
are three distinct MHD waves (fast, slow and Alfven), 
and the transfer function cannot be defined in an uni- 
vocal way. Moreover, in our simulations, unlike those of 
ICarlsson fc Steiiil (|1997l ). between the bottom boundary 
and the formation height of the Si I line there is the layer 
where sound and Alfven velocities are equal. There, the 
different modes mix up, which makes it impossible to 
find a unique transfer function. For all these reasons, 
we have decided to impose the driver at the formation 
height of the Si I line, avoiding any hypothesis about the 
propagation at deeper layers. 

Several strategies may be developed to that aim. On 
the one hand, one can set the observed oscillations as 
a boundary condition in a computational domain where 
the bottom boundary coincides with the formation height 
of the Si I. On the other hand, one can calculate the 
force which corresponds to the measured velocity and 
introduce it direct ly into the equation of motion (e.^., 
iFelipe et ani2QlQa[ ). 

In the case of the first approach, several problems arise. 
It is not valid just to set the vertical velocity, since it is 
necessary to impose at the bottom boundary the fluc- 
tuations of all the variables self-consistently. From the 
inversion of the Stokes profiles we can retrieve the varia- 
tions of all these magnitudes, but it is difficult to obtain 
reliable values with a good spatial and time resolution 
related to a si ngle layer in geometrical height, not opti- 
cal depth (see iRodriguez Hidalgo et al.li2QQ l). Another 
option is to calculate the polarization relations of all the 
variables which agree with the vertical velocity measured 
from the Doppler shift, but it is a tough work in such a 
realistic case. 

We found thus more convenient to introduce the re- 
trieved force as a source function Sz{x, zsi^t) in the mo- 
mentum equation. This driver introduces mechanical en- 
ergy in the system, and consequently, the energy equa- 
tion would also have to be modified. However, in this 
particular simulation, we use the internal energy instead 
of imposing conservation of the total energy. We have 
omitted the introduction of an additional term in the 
internal energy equation since the energy of the driver 
should not be employed in heating the plasma. The ver- 
tical force has the following dependence on the velocity: 



(3) 

where Vz,obs{x^ zsi) is the observed with the Si I A 10827 
A line velocity amplitude at frequency v for all the spa- 
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Fig. 4. — Variation of the amplitude of the velocity obtained at 
the stationary stage of simulations with the frequency of the force 
driver, normalized to the maximum value. 

tial positions x along the sht, and the functions Ai^v) and 
A(i/) include the dependence of the force amplitude and 
phase delay on frequency. The form of A(y) was evalu- 
ated by analyzing numerically the response of the atmo- 
sphere to a force driver with a fixed frequency. We have 
carried out a set of simulations using a harmonic driver 
of fixed V located at the photosphere with the force am- 
plitude independent of v. The amplitude of the velocity 
retrieved after reaching the stationary regime of the sim- 
ulations is chosen as the response of the atmosphere to a 
harmonic wave. Figure 3] shows the values obtained for 
all the simulations performed, normalized to the maxi- 
mum amplitude. The amplitude Aiv) increases almost 
linearly from v = 1.^ mHz to around 5.2 mHz, and for 
higher frequencies it decreases as the inverse of the fre- 
quency. Its maximum is located at a frequency close to 
the local cutoff frequency. The particular form of this 
dependence varies with the parameters of the sunspot 
atmosphere, and it critically depends on the location of 
the driver in height. We found that a change in the 
height location of the driver shifts the frequency of the 
maximum response. At those frequencies for which the 
response of the velocity to the force is very efficient (for 
example, at 5.3 mHz) the factor which multiplies the 
velocity to obtain the force from Equation (j3j) must be 
lower than for those frequencies with a poorer response 
(for example, at 3 mHz). For this reason A{uj) appears 
in the denominator in Equation ([3]). 

To get an idea about the form of A(i^) (phase re- 
sponse), we have carried out a simulation using Equation 
(jSj) to introduce the observed driver in the simulations, 
including the expression of A(v) from Figure H] but ig- 
noring the phase dependence A(i/). The phase difference 
between the observed velocity and the one used in the 
simulation at the height where the driver was introduced 
shows a phase delay of = 7r/2 for frequencies be- 
low 4 mHz and A(/) = for frequencies higher than 6 
mHz, with an almost linear variation between 4 and 6 
mHz. The function A(cj) has been constructed in order 
to characterize this behavior. As the velocity at low fre- 
quencies needs a quarter of a period to account for the 
variations of the force, at these frequencies A(z/) shifts 
the source driver backwards. We have set = 1.5 mHz 



and vi = 20 mHz. 

5. ANALYSIS OF THE SIMULATIONS 
5.1. Set up 

We have introduced the driver from Equation (|3]) as 
a force perturbation in the MHS model of the observed 
sunspot obtained in Section [3l The height where the 
optical depth at 5000 A is unity at the quiet Sun at- 
mosphere was chosen as z = {) Mm. According to 
iBard Carli^ (j^08^, in a model of sunspot the Si I 
line forms at a geometrical height of 308 km above the 
height where continuum rsooo = 1, which corresponds 
to z = —42 km for the adopted Wilson depression of 
350 km where the zero-height level is defined by the non- 
magnetic atmosphere. We have to take into account 
that we need at least 0.9 — 1 Mm of atmosphere above 
the location of the driver to reproduce the travel of the 
wave from the Si I line to the He I line, according to the 
geometrical height d ifference between these two layers 
(jCenteno et "all l2 0061 : ll^lipe et al.ll2010b[ ). At those high 
layers, the Alfven speed of the sunspot MHS model takes 
very high values, which requires an extremely small time 
step in the simulations. In order to save computational 
time and avoid problems with the top PML boundary, we 
have located the driver slightly deeper, at z = —100 km. 
At this height the ratio c\/v\ is around 0.8. Since we 
are introducing a vertical force in a magnetically domi- 
nating region with an almost vertical magnetic field, it 
mainly drives longitudinal waves which correspond to a 
slow magneto- acoustic mode. 

In order to compare the numerical simulation with the 
observational data, we have assigned a fixed z to the 
formation height of each spectral line, and we have as- 
sumed that the vertical velocity at that location corre- 
sponds to the velocity measured from the Doppler shift of 
the line. Fol l owing the formation heights retrieved from 
iFelipe et al.l ()2010bl ) and taking into account that the ve- 
locity obtained from the observations of the Si I line was 
imposed as a driver at zsi = —100 km, the layers of the 
computational domain selected as representative of the 
heights of formation of the rest of the spectral lines are 
zpe = 175 km, zca = 600 km and ZHe = 725 km for the 
Fel A 3969.3 line. Call H core, and He I line, respectively. 

In the horizontal directions the computational domain 
cover 14.8 x 8.4 Mm, with a spatial steps of Ax = Ay = 
100 km. In the vertical direction it spans from z = —0.6 
Mm to z = 1 Mm, excluding the PML layer, with a spa- 
tial step of Az = 25 km. We set the force Sz{i^^ zsi^t) 
for all the heights inside a layer of a chosen thickness in 
the y and z directions, but smoothly modulated to zero 
after a few grid points, and covering the slit of the obser- 
vations in the x direction. In this manner, the driver only 
acts in a narrow layer around the formation height of the 
silicon line for an elongated region in the x direction. 

5.2. Oscillations at the formation heights of the spectral 

lines 

Figure [5] shows a comparison between the LOS veloc- 
ity map of the umbra observed with the Si I hue (top 
panel) and the vertical velocity of the simulation at the 
height where the driver was introduced (middle panel). 
Negative velocities (appearing as black shaded regions) 
indicate upfiows, where the matter moves toward the ob- 
server, while white regions are downflows. Both wave 
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Fig. 5. — Velocity maps of the Sii line (the vertical direction is 
the direction along the slit, the horizontal direction is time). Top: 
Observational, measured from the Doppler shift of the Si I line; 
middle: numerical, vertical velocity at the formation height of the 
Si I line; bottom: comparison of the observational (red line) and 
numerical (blue line) velocity at x=l.l Mm. 

patterns are almost identical, with an amplitude below 
0.3 km s~^, indicating that the evaluation of the param- 
eter A{h') and A(i/) for Equation (j3]) is correct. The 
bottom panel shows the time evolution of the velocity 
at a certain location inside the umbra, confirming once 
again that the simulation fits well the observations. 

To perform a more detailed comparison between the 
velocities at the height of the driver, we show in Figure 
[6] a spectral analysis of the simulation and the obser- 
vation. The top panel illustrates the power spectra of 
the observed velocity (red dashed line) and the numer- 
ical vertical velocity at the height of the driver (blue 
solid line) averaged over slit. The ratio between the am- 
plitudes of both velocities (simulated/observed) is given 
in the middle panel. In the whole frequency range the 
ratio is around unity, indicating a good match between 
the driver and the real oscillation. The power peak at 
3 mHz is a bit lower in the simulation. For frequencies 
above 7 mHz, where the power is very low, the ratio de- 
parts slightly from unity, varying between 0.7 and 1.3. 
The bottom panel of Figure [6] shows the phase difference 
between the measured and simulated velocities. At each 
frequency we have calculated histograms of the phase dif- 
ference in all the spatial points inside the umbra. The 
color scale of the panel indicates the relative occurrence 
of a given phase shift, spanning from black (low) to red 
(high). Negative phase difference means that the sim- 
ulated velocity lags the observed one. Both oscillatory 
signals are in phase for all frequencies. 

Now we can compare the observed velocities obtained 
from different spectral lines with the simulated ones at 
heights defined in Section ISTTl 

The Fel line is formed about 280 km above the Si I 
Hue, where the driver was introduced. Figure shows the 
velocity map observed with the Fe I A 3969.3 line and the 
simulated one at the corresponding height. The Doppler 
velocity map retrieved from the Fel line is quite noisy, 
especially in this region of th e umbra where t he spectral 
line intensity is very low (^Felipe et al.ll201Qbl ). However, 
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Fig. 6. — Top: Power spectra of the observed Si I velocity (red 
dashed line) and the simulated velocity at the height where the 
driver is introduced (blue solid line); middle: ratio of the simulated 
amplitude to the observed one; bottom: differences in phase. 
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Fig. 7. Velocity maps of the Fei A 3969.3 hue. Top: Obser- 
vational, measured from the Doppler shift of the Fe I line; middle: 
numerical, vertical velocity at the formation height of the Fe I line; 
bottom: comparison of the observational (red line) and numerical 
(blue line) velocity at x=l.l Mm. 

the strongest wave fronts, with an amplitude of almost 
0.5 km s~^, can be recognized and compared with the 
simulated wave pattern, showing a good agreement. For 
example, the wavefronts around t = 20 min or the ones 
between t = 47 and t = 55 min can be clearly identified in 
the simulation, with a similar amplitude (bottom panel 
of Figure [7|). 

The power spectrum of the simulation at the forma- 
tion height of the Fel A 3969.3 line is given in Figure [8l 
The numerical simulation reproduces very well the power 
peak at 3 mHz, while the frequency of the secondary 
power peak is shifted from 5.4 mHz in the observations to 
6.2 mHz in the simulations. At frequencies above 8 mHz 
the observational power spectra shows higher power, but 
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Fig. 8. — Power spectra of the observed Fei A 3969.3 velocity 
(red dashed line) and the simulated velocity at its corresponding 
height (blue solid line), averaged over the umbra. 



we suspect that this high frequency power corresponds 
to the noise present in the velocity signal of the Fe I lines. 

Figure [9] illustrates the observational and numerical ve- 
locity maps in the case of the LOS velocity measured with 
the Call H core. The Call H core is formed at the chro- 
mosphere, and waves have propagated upwards about 
700 km from the formation height of the Si I line in order 
to reach this layer. Note that in the simulated velocity 
map (middle panel) the velocity signal is almost null dur- 
ing the first 2 minutes, due to the time spent by the slow 
magneto-acoustic waves to cover the distance between 
the driver and this height travelling at the sound speed. 
During this travel the period of the waves is reduced to 
around 3 minutes and their amplitude increases, reach- 
ing peak-to-peak values of almost 8 km s~^. The bottom 
panel of Figure [9] shows that the oscillations develop into 
shocks. This behavior is well reproduced by the numer- 
ical simulation. The simulated velocity map reproduces 
reasonably well the observed oscillatory pattern. Only in 
the temporal lapse between t = 27 and t = 40 min the 
simulated pattern differs significantly from the observa- 
tions. Most of the observed wave fronts can be identified 
in the simulation, although their spatial coverage of the 
umbra can be slightly different. For example, in the ob- 
servations at t = 50 min a wavefront shaded in white cov- 
ers from the limit of the plotted velocity map at x = 2 
Mm to around x = Mm, while in the simulations it 
extends from x = 2 Mm to almost x = 1 Mm. These 
differences may be due to the limitations in the config- 
uration of the numerical simulation: on the one hand, 
the MHS atmosphere is an axisymmetric model, which 
is obviously not the case of the real sunspot. Thus, the 
distance travelled by the waves along the field lines may 
be different, producing a phase lag. On the other hand, 
we have introduced the driver only in a region of the 
umbra along the slit of the observations, and we have 
ignored the driving of waves in the rest of the (non ob- 
served) umbra. 

During the first 20 minute of the simulations there is 
some phase shift with respect to the observations, and 
the amplitude of the simulations is lower. Note that this 
time lag is evident in the wavefronts with the largest 
amplitude, where the nonlinear it ies are clear. We will 
discuss this behavior in Section [71 

The spectral line formed the highest observed is the 
He I line, which is formed around 100 km above the 
Ca II H core. The comparison between the Doppler veloc- 
ity of this line and the vertical velocity of the simulation 
at the corresponding height is given in Figure [TOl Similar 
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Fig. 9. — Velocity maps of the Can H core. Top: Observational, 
measured from the Doppler shift of the Call H line; middle: nu- 
merical, vertical velocity at the formation height of the Call H line; 
bottom: comparison of the observational (red line) and numerical 
(blue line) velocity at x=0.9 Mm. 
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Fig. 10. — Velocity maps of the Hei line. Top: Observational, 
measured from the Doppler shift of the He I line; middle: numerical, 
vertical velocity at the formation height of the He I line; bottom: 
comparison of the observational (red line) and numerical (blue line) 
velocity at x=0.9 Mm. 



to the case of the Call H, most of the wavefronts of the 
observations can be clearly identified in the simulations, 
except in the temporal range between t = 21 and t = 40 
min. The match between the observed and simulated 
amplitudes is also remarkable. Both maps seem to be 
almost in phase. There is some phase delay which co- 
incides with the strongest shocks, but it is smaller than 
the one obtained for the Call H core. 

Only those waves with frequency above the cutoff can 
reach the chromosphere. The increase of the amplitude of 
these waves with height is larger than that of the evanes- 
cent low frequency waves, and the power spectra at the 
chromosphere is dominated by the peak at 6 mHz. For 
example, in the case of the power spectra of the He I line 
(Figure [11]) , both the observations and simulation have 
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Fig. 11. — Power spectra of the observed He i velocity (red dashed 
hne) and the simulated velocity at its corresponding height (blue 
solid line). 

their power concentrated around this frequency. The ob- 
servational power has three power peaks in the 3 minute 
band, located at 5.5, 6 and 7 mHz. The simulated power 
is concentrated at a single peak between the two high- 
est peaks of the observations. The simulated peak at 
5.5 mHz is lower than the observed one. The simula- 
tions also reproduce the power peaks at 7.7 mHz and 9 
mHz, and the low power at frequencies below 5 mHz. At 
frequencies above 13 mHz the simulated power is larger 
than the observational one. 

5.3. Propagation from the photosphere to the 
chromosphere 

The propagation properties of waves are analysed by 
means of the phase difference (Acj)) and amplification 
spectra. The value of Acj) gives the time delay between 
the oscillatory velocity signals from two spectral lines. 
We assume that the difference between them is mainly 
due to the difference of the formation height of the two 
lines. The amplification spectra are calculated as the 
ratio between the wave power at two layers. They give 
us information about variations of oscillation amplitude 
with height. The estimation of the statistical validity 
of the amplification and phase difference spectra can be 
done by calculation of the coherence spectra. This cal- 
culation s for the obs e rved d ataset were presented previ- 
ously in iFelipe et all (|2QlQb[ ) . 

In order to reproduce correctly the phase and am- 
plification spectra, the introduction of the term Qrad 
in the energy equation for the simulations (see Section 
[2]), to take into account the energy losses, has proved 
to be fundamental. This is so due to several reasons. 
Firstly, the radiative losses produce some damping of the 
waves, reducing their amplitude. Secondly, the cutoff fre- 
quency of the waves is also affected by radiative transfer 
fRoberts 1983; Khomenko et al. 2008). When the radia- 
tive timescale tr is small enough, the cutoff frequency 
is expected to decrease, compared to the adiabatic case. 
Finally, the inclusion of radiative losses also produces an 
increase of the phase difference compared to adiabatic 
case (see Figure 7 from lCenteno et aT]|2006f ). 

Figure [12] shows the phase difference and amplification 
spectra between the photospheric Si I line and the chro- 
mospheric He I line, for both observational and simulated 
velocities. From 1 to 7 mHz, where the coherence of the 
observations is high, the simulated phase difference pre- 
cisely matches the observed one, with a null phase differ- 
ence for frequencies below 4 mHz and an almost linear 
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Top: Phase difference spectra between Si I and He I ve- 



locities in the umbra. The color code shows the relative occurence 
of a given phase shift in the observations. The green crosses are 
the results of the simulation for all the spatial points. Bottom: 
Amplification spectra for the observation (black solid line) and the 
simulations (black dashed line). 

increase between 4 and 7 mHz. At higher frequencies 
the coherence of the observed phase difference is lower, 
but the simulated one keeps its linear increase. With re- 
gards to the amplification spectra, for frequencies above 
1.5 mHz the simulated spectra reproduces properly the 
observed one. The smallest frequencies show a high nu- 
merical amplification, which is possibly due to the dif- 
ficulties of the PML to damp these long period waves. 
The thickness of the PML layer should be proportional 
to the wavelength of the wave that must be absorbed, 
and the employed PML is obviously not optimized for 
such long period waves. However, since the power at 
these low frequencies is small, their overall influence on 
the simulations is negligible. 

A similar result is found between the velocity obtained 
with the Si I and the Call H core (Figure [13]), since the 
latter is formed just around 100 km below the He I line. 
The simulated phase difference is zero for frequencies be- 
low 4 mHz, and it increases at higher frequencies. It 
matches the observed phase shift between 2.5 and 9 mHz. 
At higher frequencies the observed phase difference is 
noisier, while in the frequency range between 1 and 2.5 
mHz it takes a value of tt or — tt. In the analysis of the 
observations we expressed our doubts about the reliabil- 
ity of this phase shift, so we are not surprised that the 
simulation does not reproduce it (see Felipe et al. 20 1(^. 

Figure [TH shows the phase and amplification spectra 
between two photospheric lines, the Si I and the Fel A 
3969.3 lines. The observed phase difference has high co- 
her ence between 2 and 8 mHz (Figure 12 in lFelipe et al.l 
l2QlQb). and in this frequency range the simulated phase 
delay fits the observational one reasonably well, showing 
a = for frequencies below 4 mHz and an increase 
for the higher frequencies. For frequencies below 2 mHz 
and above 8 mHz the observed phase difference spreads 
out and has lower coherence. The behavior of the sim- 
ulated amplification spectra is similar to the observed 
one between and 8 mHz, but the numerical amplifi- 
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Fig. 13. — Top: Phase difference spectra between Sii and Can H 
velocities in the umbra. Bot tom : Amphfication spectra. The for- 
mat is the same as in Figure fT2] 
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Fig. 14. — Top: Phase difference spectra between Sii and Fei 
A 3969.3 velocities in the umbra. B otto m: Amplification spectra. 
The format is the same as in Figure [12] 

cation is larger at the peak around 6.5 mHz. Note also 
that around 3.5 mHz the observed amplification shows a 
peak, not reproduced in simulations. This was previously 
seen in the power spectra in Figure [H At higher frequen- 
cies, the observed amplification has some high peaks, but 
they are not trustable due to the poor quality of the Fe I 
velocity map. 

The spectra between Fel A 3969.3 and He I lines are 
shown in Figure [15] and are similar to those in Figures [12] 
and [131 since all of them correspond to the pairs of lines 
including one photospheric and one chromospheric line. 
The numerical phase and amplification spectra match 
the observational ones in the frequency range between 2 
and 8.5 mHz. Out of this range the observational phase 
spectra are very noisy and with a lower coherence, mean- 
ing that these phase shifts are not reliable. At frequen- 
cies above 8 mHz the observational amplification seems 
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Fig. 15. — Top: Phase difference spectra between Fei A 3969.3 
and He I velocities in the umbra. Bottom: Amplification spectra. 
The format is the same as Figure [12] 

to decrease, but this tendency is not reproduced by the 
simulated one. 

6. ENERGY BALANCE 

The driver is located just above the layer where the 
Alfven speed is equal to the sound speed. At this height 
the cl/v\ parameter is 0.8. Since our driver is a vertical 
force, it mainly generates oscillations in the vertical ve- 
locity. Inside the umbral photosphere the magnetic field 
is almost vertical, and most of the energy introduced by 
the vertical force has an acoustic nature, which corre- 
sponds to a slow mode in this low-/3 region. Because of 
the vertical thickness of the driver, some part also acts 
at cs = va or below, and it produces some fast magneto- 
acoustic waves in the region just below the layer cs = va- 
These waves are partially transformed into a fast mag- 
netic mode above the height where cs = va, which is 
reflected towards deeper layers due to the gradients of 
the Alfven speed. The energy of this mode is very low. 
However, it must be taken into account that in this sim- 
ulation we are introducing the observed velocity at the 
formation height of the Si I line and, thus, the wave pat- 
tern below this layer is not reliable. In the real sunspot, 
waves propagate from deeper layers and in their upward 
propagation they reach the layer where cs = va- Some 
significative part of the energy of these waves is trans- 
formed into fast magnetic waves at this height, and the 
contribution of the fast modes in the low-/3 region must 
be larger than the one estimated in this simulation. 

The slow acoustic mode generated directly by the 
driver propagates upwards along the field hues. Accord- 
ing to Figure [6l at z = —100 km most of its power is 
concentrated in the 5 minute band, in a frequency band 
between 3 and 4 mHz. At this layer the cutoff frequency 
is Uc = 4.7 mHz, and it increases with height until it 
reaches the maximum value z^c = 6 mHz ai z = 200 km. 
It means that the oscillations in the 5 minute band in- 
troduced by the driver cannot propagate upwards, and 
they form evanescent waves which do not supply energy 
to the higher layers. Only those waves with frequencies 
above the cutoff can propagate to the chromosphere. 
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During the travel of high frequency waves to up- 
per layers their amplitude increases due to the drop 
of the density. At z = 450 km the waves de- 
velop into shocks and they reach peak-to-peak ampli- 
tudes around 8 km s~^ at the formation height of the 
Call H core {z = 600 km) and the He I line {z = 
725 km). These waves supply energy to the chromo- 
sphere, but, is this acoustic energy enough to balance 
the rad i ative losses of the chromosphere? According to 
lAvrettl (|1981f ). the radiative energy losses in the um- 
bral c hromosphere amount to 2.6 x 10^ erg cm~^ s~^, 
while iKneer et al.l (|1981[ ) estimate a similar value of 
1 — 2 X 10^ erg cm~^ s~^. These losses should be bal- 
anced by some energy source. Figure [16] illustrates the 
variation of the acous tic flux with height, w hich was cal- 
culated following Bra v~fc LoughheadI ( 19741 ): 



(Pivi), 



(4) 



where pi and vi are the fluctuations in gas pressure and 
velocity, respectively. It was averaged in time for all 
points inside the umbra at the stationary stage of simula- 
tions. At the position of the driver (z=-100 km) it is zero 
because it represents the average between the upward 
and downward fluxes. It grows up to 10^ erg cm~^ s~^ 
just above the driving layer and decreases with height 
due to the dissipation and radiative energy losses. The 
decrease between z = 200 km and z = 600 km is mainly 
due to the impossibility of the waves with frequencies 
below the cutoff frequency to propagate to higher lay- 
ers. Above z = 600 km the decrease of the acoustic flux 
is even more pronounced because of the dissipation pro- 
duced by the shocks and the losses caused by the low 
radiative relaxation time at chromospheric heights (we 
set = 10 s above 600 km). 

The magnetic flux associated to the fast 
magneto-acoustic waves can be calculated following 
iBrav fc LoughheadI (jl974D as 



mag 



(Bi X (vi X Bo)//io) 



(5) 



It shows a maximum of 10^ erg cm~^ s~^ at z 
km, and it drops at upper layers with a value of 
2 X 10^ erg cm~^ s~^ at z = 750 km since these waves 
are reflected down towards the photosphere due to the 
gradients of the Alfven speed (iRosenthal et al.l 120021 : 
iKhomenko fc ColladosI [2QQ6I : iFelipe et al.ll2010aD . How- 
ever, it must be taken into account that most of the 
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Fig. 17. — Acoustic energy flux map at ZHe- Positive values 
(white) indicate upward flux. 

energy of the vertical force that we introduce as a driver 
mainly goes to the slow magneto-acoustic mode, and this 
calculation of the flux carried by the fast mode may be 
underestimated. 

The acoustic flux shows important variations with 
time. For example. Figure [TTI shows the temporal evolu- 
tion of the acoustic flux inside the umbra at the formation 
height of the He I line. Positive values (shaded in white) 
represent upward flux. The wavefronts of the slow acous- 
tic waves in the 6 mHz band are accompanied by acous- 
tic energy which reaches almost 5 x 10^ erg cm~^ s~^ for 
the strongest shocks. The average value of the acoustic 
flux at this height is 3 x 10^ erg cm~^ s~^, which is an 
order of magnitude below the required value to balance 
the radiative losses. Only at the moments when shocks 
reach the chromosphere the energy supplied by the slow 
acoustic waves is similar to the chromospheric radiative 
losses. 

7. DISCUSSION AND CONCLUSIONS 

We have presented numerical simulations which closely 
reproduce the real wave propagation observed in the um- 
bra of a sunspot from the photosphere to the chromo- 
sphere. The anal ysis of these observa tions was previ- 
ously report e d in iFelipe et al.l (|2010b[ ). As proved in 
IFelipe et al.l (j2010aD . our code is able to manage the 
propagation of waves with long realistic periods, even in 
the hard conditions imposed by the high Alfven speeds 
at the sunspot chromosphere. 

The similarity between the simulations and the ob- 
servations is achieved in two steps. Firstly, we had to 
construct a MHS model of a sunspot with properties as 
close as possib le to the observed one. We ha ve followed 
the method of IKhomenko fc ColladosI (|2008l ). introduc- 
ing some modiflcations in order to take into account 
thermodynamic and magnetic properties of the partic- 
ular observed sunspot, as retrieved from the inversion of 
the Si I line. Secondly, this static model was perturbed 
by a force which drives a vertical velocity similar to the 
LOS velocity measured from the Doppler shift of the Si I 
line. 

Once the observed velocity is introduced in the MHS 
atmosphere, it mainly drives slow magneto- acoustic 
waves in the low-/3 region. Most of their power is concen- 
trated in the 5 minute band, with frequencies between 3 
and 4 mHz, and they form standing waves due to the 
higher cutoff frequency of the atmosphere. The driver 
also generates high frequency waves, which propagate 
upwards to the chromosphere. In their travel through 
the sunspot atmosphere, they reach the formation height 
of several observed spectral lines. The simulated velocity 
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maps and power spectra at their corresponding formation 
heights reproduce reasonably weh those obtained for the 
chromospheric Call H core and He I hne. In the case of 
the Fel A 3969.3 hne, the comparison is hindered due to 
the poor quahty of the observed velocity map. 

The comparison between the simulated and observed 
velocity maps at chromospheric heights reveals some 
phase delay at the strongest wave fronts, in the sense that 
the observed wave fronts lag the simulated ones. In the 
case of the Call H core the phase lag is OAtt. The delay 
of the observed relative to the simulated He I velocities is 
somewhat smaller. This delay, together with the larger 
amplitude of the observed strongest shocks, may indi- 
cate that in the dynamical sunspot chromosphere, the 
nonlinear waves shift the formation height of spectral 
lines to higher layers. This is especially evident for the 
Call H core. In the analysis of these simulations, we have 
chosen a certain height representative of the layer from 
which the information about the velocity measured with 
the observed line comes. Obviously, this procedure is a 
simplification, since the contribution function from every 
spectral line spans over a thick layer. The height that we 
ar e considering corresp onds to the average height derived 
in iFelipe et al.l ()2QlQbl ) . This approximation seems to be 
good for the photospheric lines, but not for the highly 
nonlinear chromospheric region. 

Our simulations reproduce the phase and amplification 
spectra between several pairs of lines with a remarkable 
match. For those spectra between a photospheric and 
a chromospheric signal, the phase difference shows sta- 
tionary waves with Acj) = below 4 mHz. At higher 
frequencies the waves progressively propagate with Acj) 
increasing with the frequency. The empirical cutoff of 4 
mHz is assigned from the estimation of the frequency at 
which the phase difference starts to depart from Acf) = 0. 
In the solar atmosphere, the cutoff is a local parameter 
which depends on the temperature and is stratificated 
with height. In our MHS model its maximum value is 6 
mHz, and only waves with frequencies above this value 
can propagate all over the atmosphere. However, in most 
of the atmosphere the local cutoff is below 6 mHz, and 
for this reason the phase difference spectra show propa- 
gation for frequencies in the range between 4 and 6 mHz. 

The only strong discrepancy between observations and 
simulations is found in the amplification spectra at fre- 
quencies above 8 mHz (Figures [13] and [Hj). The am- 
plification is cleary higher in the simulations than in the 
observations, and the observed power spectra of the chro- 
mospheric lines show lower power than the simulated 
equivalent (Figure [TT]) . There are several causes that 
could be responsible for this discrepancy. The compar- 
ison of the power spectra of the velocity retrieved from 
the Si I line and the one obtained in the simulation at its 
formation height (middle panel of Figure [6]) shows some 
slight differences which are especially significant above 8 
mHz. Since the amplitude of the waves increases expo- 
nentially with height, these small discrepancies could be 
the origin of the larger ones that appear at larger heights. 
For example. Figure [6] shows that the simulation has an 
excess of power at the height of the driver at 13 mHz 
and 17 mHz, and these frequencies also present higher 
power than the observations at the chromosphere (Fig- 
ure [TT]). However, some frequencies where the power of 
the photospheric velocity in the simulation is not larger 



than the power of the Si I line also present a significative 
difference at higher layers. It is possible that most of the 
power generated by the driver at high frequencies (which 
match the power obtained in the observations) does not 
correspond to real photospheric oscillations, but it is in- 
troduced by some observational limitations. In this case, 
the simulation would propagate upward some oscillatory 
power which is not present in the Sun, and although 
it vanishes at deep layers, its contribution in the chro- 
mosphere is important. Another possible cause of the 
different power of the high frequency waves at the chro- 
mosphere is the simple energy exchange implemented in 
this simulation. A more realistic treatment could pro- 
duce a stronger radiative damping of the high frequency 
waves. 

The analysis of the energy flux carried by magnetoa- 
coustic waves to high layers reveals that it is not enough 
to heat the chromosphere. The average energy supplied 
by slow acoustic waves is around 3 to 9 times lower than 
the amount required to balance the radiative losses, de- 
pending on the reference taken for the umbral chromo- 
spheric radiative losses. However, the simulation over- 
estimates the power of the waves with frequencies be- 
tween 8 and 20 mHz at the chromosphere, and, thus, 
the contribution of waves to chromospheric heating in 
sunspot could be lower. In our simulations we have re- 
jected the wav es with frequencies above 20 mHz. Ac- 
cordin g to Foss um fc Ca rlsson (2005) and lCarlsson et all 
(|2007l ) the power of waves with 20 mHz frequency can 
be assumed to be negs:ligs:ible for chrom ospheric h e ating , 
but lWedemever-Bohm etal] (j2007[ ) and iKalkofen] (j2007n 
have pointed out that this question depends on the spa- 
tial resolution of the data used. 

Previous works have found that the acoustic wave en- 
ergy is too low, by a factor of at least ten, to balance 
the rad iative losses in the non-magnetic solar chromo - 
sphere (|Fossum Carlsson] [2005] : iCarlsson et"al] 120071 ) . 
In the context of small-sc ale magnetic flux concentra- 
tions, IVigeesh et al.l ()2009[ ) obtained that the acoustic 
energy flux of transversely, impulsively excited waves is 
insufficient to balance the chromospheric radiative losses 
in the network. The fraction of the required energy sup- 
plied by acoustic waves in the magnetized atmosphere of 
a sunspot seems to be larger than in the quiet Sun. The 
chromospheric energy losses in umbral regions are 3-5 
times lower than in quiet Sun, but the energy supplied 
by the waves seems to be similar. At the photosphere, the 
surface amplitude of solar p-modes in magnetic regions 
are reduced below t hose in magnetically q u iet regions in 
the 5 minute band (|Woods fc Cramlll98l[ ). iBrown et al.l 
(|1992[ ) showed that the supression is frequency depen- 
dent, peaking at 4 mHz. The reduction of the amplitude 
occurs for frequencies below 5 mHz, while the ampli- 
tude of waves above 5.5 mHz is enhanced. Since the 
waves that propagate energy to higher layers are those 
above the cutoff frequency, the reduction of p-modes in 
regions of magnetic activity does not produce a lower ef- 
ficiency of acoustic energy propagation in sunspots. Sev- 
eral works have pointed out that the fast magnetic mode 
in the low-/3 region is reflected back toward the photo- 
sph ere and it is unable to supply energy to higher lay- 
ers ([Rosenthal et al.l 120021 : iKhomenkofc Colladosl 120061 : 
IFelipe et al.ll2010a] ). According to lCallvl ('2005'), the fast- 
to-fast conversion is more efficient for higher frequencies 
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and, thus, waves with frequencies above the cutoff are 
significantly transformed into fast magnetic waves above 
the height where va = cs and cannot carry energy to 
the chromosphere. In our simulation the energy fiux 
of the fast magneto-acoustic mode at the chromosphere 
amounts to 2 x 10^ erg cm~^ s~^, being 2-3 orders of 
magnitude lower than the acoustic fiux. The LOS ve- 
locity oscillations measured with the Si I line correspond 
to the slow acoustic waves in the low-/3 region, after the 
transformation of the upward propagating waves. Their 
power has already been reduced by the transformation to 
the fast magnetic mode. Figure [iBl shows that at the for- 
mation height of the Si I line, the average acoustic energy 
fiux is around 10^ erg cm~^ s~^, so at this photospheric 
height the energy contained in the form of acoustic fiux 
is of the order of magnitude of the amount required by 
the chromospheric radiative losses. 

As a further step in our work, we aim to synthetize the 
spectral lines from the numerical simulation in order to 
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develop a better comparison with the spectropolarimetric 
observations. This process is relatively simple in the case 
of the photospheric lines, since they are formed at LTE, 
but the non-LTE formation of the chromospheric lines 
hinders their synthesis. The last case is especially inter- 
esting to clarify the issue of the variation in the velocity 
response height of the atmosphere to nonlinear perturba- 
tions. On the other hand, we also plan to perform similar 
simulations using a bidimensional velocity field obtained 
with IBIS or HMI/SDO as a driver. 
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